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o-RuCla has been proposed recently as an excellent playground for exploring Kitaev physics on a two- 
dimensional (2D) honeycomb lattice. However, structural clarification of the compound has not been completed, 
which is crucial in understanding the physics of this system. Here, using ab-initio electronic structure calcula¬ 
tions, we study a full three dimensional (3D) structure of o-RuCH including the effects of spin-orbit coupling 
(SOC) and electronic correlations. Three major results are as follows; i) SOC suppresses dimerization of Ru 
atoms, which exists in other Ru compounds such as isostructural Li 2 Ru 03 , and making the honeycomb closer 
to an ideal one. ii) The nearest-neighbor Kitaev exchange interaction between the jeff=1/2 pseudospin depends 
strongly on the Ru-Ru distance and the Cl position, originating from the nature of the edge-sharing geometry, 
iii) The optimized 3D structure without electronic correlations has P31m space group symmetry independent 
of SOC, but including electronic correlation changes the optimized 3D structure to either C2 /?ti or Cmc2i 
within 0.1 meV per formula unit (f.u.) energy difference. The reported P3i 12 structure is also close in energy. 
The interlayer spin exchange coupling is a few percent of in-plane spin exchange terms, confirming a-RuCla is 
close to a 2D system. We further suggest how to increase the Kitaev term via tensile strain, which sheds new 
light in realizing Kitaev spin liquid phase in this system. 


I. INTRODUCTION 

There have been a number of studies on quasi-two- 
dimensional systems having both spin-orbit coupling (SOC) 
and on-site Coulomb interactions, which are believed to host 
unconventional magnetic orders and spin liquid phases 
One promising candidate is a-RuCls, where edge-sharing 
RuCle octahedra form two-dimensional RuCls layers in 
which Ru honeycomb layers reside^“^^ Compared to its 5d 
transition metal oxide counterparts a-A 2 lrOs (A=Li,Na)^^“^^, 
a-RuCls has closer-to-ideal RuCle octahedra^, so it was pro¬ 
posed as an excellent platform to explore the Kitaev physics 
and related magnetism despite weaker A few 

recent reports suggest the presence of strong Kitaev-type 
bond-dependent exchange interactions in a-RuCls^, which 
originate from the cooperation between the intermediate SOC 
in Ru atom and the Coulomb interaction^. A zigzag-type mag¬ 
netic order within the RuCls layer is also predicted and ob¬ 
served, which is proximate to the Kitaev spin-liquid phase^’^. 

In previous studies a-RuCls was considered as a two- 
dimensional system with an ideal Ru honeycomb lattice, 
but such assumption needs further clarification. A poten¬ 
tial Ru layer distortion, which is observed in an isostruc¬ 
tural compound Li 2 Ru 03 ^^’^^, might happen in this com¬ 
pound. Furthermore, a-RuCls has a three-dimensional crys¬ 
tal structure consisting of RuCls layer stacking, and inter¬ 
layer coupling and interaction terms can introduce another 
complication. Experimentally, both P3il2 and (72/m space 
groups have been suggested as the crystalline symmetry in 
this compound^’^’^^’^^’^^. As an illustrative example, Fig. 1(a) 
shows the crystal structure of a-RuCls with a (72/m space 
group symmetry, where adjacent RuCls layers within the unit 
cell is related to each other by a translation along the a-axis in 
the figure. Stacking faults can easily be introduced in this lay¬ 
ered structure as in the case of a-A 2 lrOs^^, which obscures 
further clarification of the crystal structure. Effect of inter¬ 
layer exchange interactions from the layer stacking on the 



FIG. 1. (Color online) (a) Crystal structure of o-RuClawith (72/m 
space group. Solid lines depict a monoclinic unit cell, (b) Schematic 
view of three triangular sublattices on which Ru and Cl layers are 
located. Stacking indices for Ru honeycomb and Cl triangular layers 
are shown on the right side of (a), where indices for Ru and Cl layers 
are expressed as capital and lowercase letters respectively. 


ground state magnetic properties of this system is not well 
understood either. More interestingly, a sample-dependent 
two-transition behavior is reported, where two different mag¬ 
netic order peaks at T/vi — 14 K and Tn 2 — 8 K with two- 
and three-layer c-axis periodicity, respectively, are observed in 
neutron diffraction measurement^ 7 These issues pose a ques¬ 
tion on the relation between crystal structure and magnetism 
in this system. 

In pursuit of such motivations, in this work we perform ab- 
initio calculations for the structural properties of a-RuCls and 
their impact on magnetism. We present three main results; i) 
Role of SOC and zigzag magnetic order on the single-layer 
RuCls structure is discussed. We found that SOC prefers 
ideal honeycomb lattice by preventing Ru-Ru dimer forma- 
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tions, and the presence of in-plane zigzag order tends to give 
small monoclinic distortion commensurate to the magnetic or¬ 
der. ii) Effect of Ru-Cl and Ru-Ru distance to the exchange 
interactions and magnetism is discussed, where the hopping 
channels within the nearest-neighbor(NN) Ru t 2 g orbitals and 
the resulting exchange interactions between the SOC-induced 
jfgf^=l/2 pseudospins strongly depend on the Ru-Cl and Ru- 
Ru distance. Such behavior originates from the existence of 
multiple hopping channels in the t 2 g orbitals, which enables 
'leveraging' magnetism with rather small amount of struc¬ 
tural changes, iii) Stability of crystal structures with different 
stacking orders is discussed by comparing relative total ener¬ 
gies. We have found that, structures with and Cmc2i 

space group symmetries are most favorable with almost de¬ 
generate energies. Previously suggested P3il2 structure^’^^ 
yields total energy comparable to those of C2/m and Cmc2i 
structures with the energy difference smaller than 0.4 meV per 
formula unit (f.u.). Energy differences between different in¬ 
terlayer magnetic orders are smaller than 0.1 meV / f.u., and 
magnitude of interlayer exchange interactions estimated from 
interlayer hopping integrals are smaller than 0.05 meV. These 
observations justify the employment of two-dimensional spin 
models in exploring magnetism in a-RuCls. We further pro¬ 
pose how to increase the Kitaev term using tensile strain or 
uniaxial pressure to realize the Kitaev spin liquid phase. 

This manuscript is organized as follows. After show¬ 
ing computational details in Sec. II, structural properties of 
single-layer RuCls and its relation to magnetism is presented 
in Sec. III. The effect of SOC and zigzag magnetic order to 
the single-layer RuCls structure, and the relation between the 
structure and magnetism are discussed in Sec. Ill A and Sec. 
IIIB, respectively. In Sec. IV and V, results on the stacking 
without and with the Coulomb interaction and magnetism are 
shown, respectively. Summary and conclusion follow in Sec. 
VI. 


II. COMPUTATIONAL DETAILS 

Eor the electronic structure calculations, we employed the 
Vienna ab-initio Simulation Package (VASP), which uses the 
projector-augmented wave (PAW) basis set^"^’^^. 370 eV of 
plane wave energy cutoff was used, and for /c-point sampling 
15x15 and 8 x 6 x 4(6) Monkhorst-Pack grid were adopted for 
single-layer primitive cell and monoclinic cells with three 
(two) layer c-axis peroidicity. Tetrahedron method with 
Blochl correction was used for the calculation of density of 
states^^. On-site Coulomb interactions are incorporated us¬ 
ing the Dudarev’s rotationally invariant DET-^f/ formalism^^ 
with effective f/eff = U — J = 2 eV. Eor each configura¬ 
tion with different unit cell, f/eff value, and magnetic order, 
structural optimization is performed with a force criterion of 1 
meV / A. Unless specified, a revised Perdew-Burke-Emzerhof 
generalized gradient approximation (PBEsol)^^ was used for 
structural optimization and total energy calculations. Note 
that, PBEsol functional yielded reasonable results for the 
stacking order of bilayer transition metal dichalcogenides 
in comparison to the van der Waals functionals^^. Results 


with employing vdW functionals are shown in Appendix. 
Interlayer hopping integrals were obtained by employing 
maximally-localized Wannier orbital (MLWE) formalism^^’^^ 
implemented in Wannier90 package^^. Also, for compari¬ 
son of the magnetism in the single-layer structures in Sec. 
Ill, a linear-combinaion-of-pseudo-atomic-orbital basis code 
OPENMX^^’^^ was used, where double zeta plus polarization 
(DZP) bases, 500 Ry of energy cutoff for real space integra¬ 
tions, and the Perdew-Zunger parameterization for the local 
density approximation were employed^^’^^. 


III. RELATION BETWEEN STRUCTURE AND 
MAGNETISM IN RUCL3 SINGLE LAYER 

In this section, structural changes due to the lattice opti¬ 
mization and their effect to the magnetism is discussed in the 
RuCls single layer. The initial trial structure we chose is the 
one reported in Ref. 37, which was used in the Ref. 8. The lat¬ 
tice optimization gives rise to in-plane structural changes, and 
here we present the optimized structures focusing on the dif¬ 
ference from the old one. Since we found that such behavior 
and the resulting changes in magnetism also occur in the full 
3D structures, which are presented in Sec. IV and V, below 
we first discuss the single layer results. 


A. Effect of SOC on in-plane Ru dimerization 

Eirst, the effect of SOC and magnetism with Ueff on a Ru 
honeycomb lattice is discussed in this subsection. Eig. 2 
summarizes the results, where the sizes of Ru displacements 
6 from the ideal honeycomb lattice after structural optimiza¬ 
tions under different conditions are shown. Positive and neg¬ 
ative values of 6 in Pig. 2(b) correspond to Ru dimerization 
and Ru zigzag chain formation, respectively, as shown in Pig. 
2(a). Note that, the lattice constants are fixed to the exper¬ 
imentally observed a = oq = 5.96A and b=\/3ao. Without 
including SOC and Coulomb interactions, the two Ru atoms 
in the unit cell tend to dimerize to lower the energy as shown 
in Pig. 2(a). The presence of dimer formation is robust 
against different choice of exchange-correlation functionals 
— Perdew-Zunger parametrization of local density approxi¬ 
mation (LDA)^^, PBE^^, and PBEsol — with slightly different 
size of 6 as shown in Pig. 2(b). Similar dimer formation was 
reported in other layered honeycomb compound Li 2 Ru 03 , of 
which origin is suggested to be the cr-like direct bonding be¬ 
tween the neighboring Ru t 2 g orbitals 

Since the dimer formation breaks the Ru t 2 g degeneracy, we 
expect that SOC would not favor the dimer formation. The 
spin-orbit entangled jeff orbitals, which emerges under the 
presence of cubic crystal fields and SOC"^®’'^^ does not favor 
orbital polarization between the t 2 g — dxy, and dyz — 
orbitals. Indeed, structural optimizations including SOC yield 
significant reduction of dimerization as shown in the middle 
of Pig. 2(b). Although there are small differences between 
EDA, PBE, and PBEsol results, the role of SOC in preventing 
the dimerization is evident. Additionally, inclusion of the on- 
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FIG. 2. (Color online) (a) Schematic figure of Ru honeycomb with 
colored circles depicting Ru sites. Grey dotted and black dashed 
squares represent the primitive and monoclinic unit cells respec¬ 
tively, where colors on Ru sites show the zigzag magnetic order in 
a monoclinic unit cell, (b) Size of Ru distortion 6 under differ¬ 
ent exchange-correlations functionals and with/without presence of 
SOC, t/eff, and in-plane zigzag magnetic order. Note that, positive 
and negative S correspond to Ru dimer and zigzag chain formations, 
respectively. 
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site Coulomb interaction without the presence of magnetism 
is expected to enhance the idealness of the Ru honeycomb lat¬ 
tice, since it was shown previously that the on-site Coulomb 
interaction effectively enhances the size of SOC^’"^^. 

Next we show the effect of in-plane zigzag magnetic or¬ 
der, which is predicted to occur when SOC and the Coulomb 
interaction are incorporated in ab-initio calculations^ and ob¬ 
served in experiments^’Right columns of Fig. 2(b) show 
the results from calculations including SOC, = 2eV, and 
the zigzag order. The enlarged monoclinic unit cell and the 
magnetic configuration are shown in Fig. 2(a), where the red 
and blue colored circles represent Ru sites with antiparallel 
moments to each other. Regardless the choice of functional, d 
shows negative values with almost same magnitude. The re¬ 
sulting structure is commensurate to the zigzag magnetic or¬ 
der as shown in Fig. 2(a), suggesting a finite magneto-elastic 
coupling in this compound. 


B. Effects of Cl displacement and lattice constant change to 
the exchange interactions between the jeff=1/2 pseudospins 

Here we discuss the Cl displacement after the optimization 
and its impact to the exchange interactions between the neigh¬ 
boring Ru jfeff=l/2 pseudospins. Fig. 3 shows the displace¬ 
ment of Cl atoms after structural optimization, where the two 



FIG. 3. (Color online) Schematic figure showing the direction of Cl 
displacement from the ideal position after the structural optimization. 
Three inequivalent NN bonds — Z-, X-, and Y- bonds — and the 
displacements of participating Cl atoms therein are depicted by red, 
blue, and green planes and arrows, respectively. Cl triangles located 
above and below the Ru plane are represented as solid and dotted 
triangles, respectively. 


Cl atoms participating in each NN Ru bond move toward the 
bond center. When the in-plane lattice constants are fixed to 
be a = ao and b=^/3ao, structural optimization with SOC only 
(no (7eff and magnetism) yields reduced Cl height of 1 .43 A to 
1.34 A with respect to the Ru plane, and the Cl triangles above 
and below Ru plane rotates by 2.1° in opposite direction as 
shown in the figure. The Ru-Cl-Ru NN bond angle increases 
from 89.1° to 93.8°. After allowing the lattice constants to re¬ 
lax, the lattice constants reduce to a = 0.98lao and b = 0.9866o 
when SOC was employed with the monoclinic distortion al¬ 
lowed. With (7eff = 2 eV and the zigzag magnetic order, they 
are increased to a = 1.01 lao and b = I.OO 660 . The averaged 
Ru-Cl distance changes from 2.34 A to 2.36 A in the nonmag¬ 
netic calculation with (7eff = 0 eV to the magnetic results with 
= 2 eV, but both of them are shorter than the distance of 
2.45A in the initial trial structure. Note that, when the mono¬ 
clinic distortion is allowed, the NN Z-bond in Fig. 3 becomes 
inequivalent to the X and Y bonds, where the X and Y bonds 
form the zigzag chain in Fig. 2(a). Also, no Ru-Cl bond length 
disproportionation is observed in all of our results, implying 
no Jahn-Teller distortion in this system. 

Due to the presence of inversion symmetry at the bond cen¬ 
ter and additional trigonal distortion in RuCle octahedra, the 
hopping integrals between the NN Ru t 2 g Wannier orbitals 
have the following form^^’"^^, 

fh t2 U\ 

T = t 2 ^4 , 

V u ^4 h ) 

where each hopping channel is displayed in Fig. 4 with the 
participating Ru t 2 g Wannier orbitals therein. As shown in the 
figure, while ti originates mainly from the 5- and cr-like d-d 
direct overlap integrals, ^2 is mostly from the 7r-type indirect 
overlap dominated by d-p-d hopping between the Ru and in- 
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FIG. 4. (Color online) Four major NN hopping channels — (a) ta, 
(b) t 2 , (c) ti, and (d) — within the t 2 g subspace. For each hopping 

channel, the participating t 2 g Wannier orbitals are plotted, where the 
schematics for each channel is represented in the inset. Note that, ts 
and t 2 terms depends more sensitively to the structural change than 
ti and t 2 . 


tervening Cl p orbitals. Note that, ts channel has both the d-d 
direct overlap and d-p-d indirect overlap which has opposite 
signs to each other. Also, due to the small trigonal distortion, 
the small and terms are introduced, where the difference 
between them introduced by the monoclinic distortion is neg¬ 
ligibly small. 

Table I shows the hopping terms from the Wannier orbitals 
for four crystal structures optimized with different conditions. 
There are the old P3il2 structures? used in previous work, 
structure with internal coordinates and lattice constants opti¬ 
mized with SOC, structure with only internal coordinated op¬ 
timized (fixed a=ao and b=bo), and the one optimized with 
SOC, Peff and the zigzag order. Hereafter we denote the 
structures as Case 0 to III, respectively, as stated in Table I. 
With those optimized structures, calculations of the Wannier 
orbitals were performed without the inclusion of SOC, Peff. 
and magnetism. Surprisingly, the hopping integrals are show¬ 
ing huge dependence to the structural change. Especially, the 
ts term varies from -0.229 to -0.062 eV depending on the 
structures, and t 2 also varies from 0.114 to 0.191 eV. Com¬ 
paring the Case 0 and II results, the effect of Cl relaxation is 
to enhance ^2 and suppress The effect of increasing Ru- 
Ru distance, which can be seen by comparing Case I to III, is 
also similar to the role of Cl relaxation with less dramatic but 
still substantial trend. Such tendency can be understood from 
the character of participating Wannier orbitals shown in Fig. 
4. The ts term, the most sensitive to the structural change, 
originates from the two distinct channels; one from the cr-like 
direct d-d overlap and another from d-p-d indirect channel. 
The two channels has opposite sign to each other, with minus 
sign for the d-d channel and plus sign for the d-p-d channel. 
As a result, enhancing d-p-d channel by reducing the Ru-Cl 
distance or increasing the Ru-Cl-Ru angle will lead to bet¬ 


ter cancellation of the dominant d-d channel and reduction of 
the overall ts term as shown in Table I. Enhancement of t 2 
after Cl relaxation is also easy to understand since it mostly 
comes from the 7r-like d-p-d channel, while the ts dominated 
by the (5-like d-d channel is reduced as the Ru-Ru distance is 
increased. The trend for the small ^4 term is less clear, but it 
tend to enhance when there are more trigonal and monoclinic 
distortion. 

From the NN t 2 g hopping terms, one can estimate the val¬ 
ues of exchange interaction terms in the jeff = 1/2 spin 
Hamiltonian 


{ij) 

where the bond-dependent 3x3 matrix ^ has the form of 

(j V r \ 

M = r j r' . 

\V' v J^K ) 

Note that, undergoes simultaneous cyclic permutations 
of rows and columns depending on NN bond directions. Ex¬ 
plicit expressions for the Heisenberg J, the Kitaev K, and 
the symmetric anisotropy terms T and T' in terms of the hop¬ 
ping integrals, U, and the Hund’s coupling Jh are reported 
in Ref. 16 and 43. Using the values of ti listed in Table I 
and setting U = 3eV and JhIU = 0 . 2 , we can calculate the 
values of exchange interactions which are listed in Table I. 
Note that, changing the values of U and JhIU changes does 
not change the ratio between the exchange interactions when 
Juju > 0.05. As shown in the table, among the exchange 
interactions, the Kitaev term shows dramatic change of chang¬ 
ing sign after the Cl relaxation. This is due to the enhancement 
and suppression of t 2 and ts terms. Increasing Ru-Ru distance 
gradually enhances K and reduces J and T, so driving the sys¬ 
tem closer to the Kitaev spin liquid limit with ferromagnetic 
K. Note that, comparing Case II and III, increasing the lat¬ 
tice constant by 1 % enhances the K term significantly. This 
implies the possibility of controlling the magnetism and real¬ 
izing the Kitaev spin liquid phase with rather small amount of 
structural change such as epitaxial strain or uniaxial pressure. 
Another noticeable feature is the small but non-negligible T' 
term from the trigonal distortion, which can stabilize the ex¬ 
perimentally observed zigzag order near the Kitaev spin liquid 
phase with K < 0^^ . 

Finally, we discuss the evolution of the magnetic mo¬ 
ments direction in the zigzag order with respect to structural 
changes. Fig. 5(a) shows the schematic figure of the zigzag 
order with an angle of the moments 0 with respect to the a- 
axis. Note that, in all of our calculations the moments were 
residing on the ac-plane. In the Case 0 structure, both in 
the OpenMX and VASP results, the moments were paral¬ 
lel/antiparallel to the a-axis {i.e. 0 = 0 ), consistent to our pre¬ 
vious result^. After the structural optimization the moments 
gain nonzero 0, which tends to increase when the lattice con¬ 
stant increases as shown in Fig. 5(b). There is difference in 
0 between the results from the two different codes, but the 
tendency of increasing angle remains the same. We speculate 
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javg 

“Ru-C1 

C^Ru —Ru 

(in A) 

ti t2 ts t4 

(in eV) 

j K r r' 

(in meV) 

Case 0 structure: old P3i 12 structure 
(from Ref. 3, a=ao, b=bo) 

NN 2.45 

3.44 

+0.066 +0.114 -0.229 -0.010 

-3.50 +4.60 +6.42 -0.04 

Case I structure: a=0.981ao, b=0.986bQ 
(structure optimized with SOC) 

NN-Z 

3.40 

+0.058 +0.177 -0.154 -0.022 

-2.67 -4.52 +1.21 -0.67 

2.34 

NN-X/Y 

3.38 

+0.060 +0.165 -0.160 -0.018 

-2.81 -3.07 +6.99 -0.47 

Case 11 structure: a=ao, b=bo 

(structure optimized with SOC and lattice constants fixed) 

NN-Z 

3.44 

+0.044 +0.178 -0.109 -0.019 

-1.49 -6.71 +5.28 -0.69 

2.36 

NN-X/Y 

3.44 

+0.042 +0.176 -0.107 -0.030 

-1.55 -6.47 +5.24 -1.08 

Case III structure: a=1.011ao, b=1.006bo 
(structure optimized with SOC, Ceff, and zigzag order) 

NN-Z 

3.47 

+0.036 +0.191 -0.062 -0.024 

-0.74 -9.34 +3.71 -1.04 

2.36 

NN-X/Y 

3.47 

+0.037 +0.182 -0.075 -0.026 

-1.09 -7.64 +4.38 -0.87 


TABLE L Values of the averaged Ru-Cl distances, NN Ru-Ru distances, hopping integrals, and examples of exchange interactions for U = 
3 eV and Jn/U = 0.15^^. Case I to III structures were optimized with different conditions stated inside the table, while Case 0 structure is 
from Ref. 37. In Case II, lattice constants are fixed to be ao and bo, while in Case I and III they are allowed to relax. Hopping integrals and 
exchange interacitons are shown in eV and meV units, respectively. For comparison, the values of hopping integrals and exchange interactions 
from the old structure (Case 0) in Ref. 8 are listed. 


that such behavior may originate from the Cl relaxation and 
the resulting change of exchange interactions, especially the 
change of the ratio between K and F terms. Also, as the lat¬ 
tice constant is enlarged, the two zigzag chains with antiparal¬ 
lel moments in the unit cell begin to develop the difference in 
0, so resulting net ferromagnetic component in the ac-plane. 
The origin of such behavior is unclear at this point. 


IV. STACKING WITHOUT Ues AND MAGNETISM 

Next let us study the stacking order of RuCls. First we dis¬ 
cuss their relative total energies without including t/eff and 
magnetism. As mentioned in Sec. II, here we show the re¬ 
sults with using PBEsol functional, and their comparison with 
vdW functional calculations are shown in the Appendix. Note 
that, PBEsol results give the same lowest energy configura¬ 
tions with other vdW results, and the closest c-axis constant 
to the experimentally observed one as welP. 

Fig. 6 shows five unit cells we considered in this work, 
where the upper and lower panels show the side view of unit 
cells and top view of Ru honeycomb layers respectively. Once 
we consider Ru honeycomb as a triangular layer by ignoring 
Ru hollow sites, a-RuCls crystal structure can be understood 
as a stacking of Ru and Cl triangular layers with three triangu¬ 
lar sublattices (a/A, b/B, and c/C, where capital and lowercase 
letters denote Ru and Cl layers respectively) shown in Fig. 
1(a) as a degree of freedom. In Fig. 6, each different structure 
can be understood as a sequence of sublattice indices. Note 
that, within a RuCls layer, any two Ru or Cl layers cannot be 


in a same sublattice. As we take into account Ru hollow sites, 
additional degree of freedom is introduced to each Ru layer, 
and we denote this with primes in the triangular sublattice in¬ 
dex (for example. A, A’, and A” as shown in the figure). 

For structures with three-layer c-axis periodicity, we choose 
unit cells with P3il2 and C2/m space groups. Note that, the 
C2/m structure was reported also as the space group of this 
compound^’^^, and is similar to the P3il2 structure. The ma¬ 
jor difference in two structures is the c-axis ordering of the 
Ru honeycomb layers, where in the C2/m unit cell three Ru 
layers are related by translation by (a + a + c)/3 while in 
the P3il2 cell they are related by threefold screw axis. Be¬ 
sides, since the neutron diffraction result identified a magnetic 
peak with two-layer c-axis periodicity atT/vi = 14Kina 
poly crystalline sample^ \ we consider two-layered unit cells 
as well. Avoiding two Cl“ triangular layers belonging to adja¬ 
cent RuCls layers to locate on top of each other (i.e. sitting on 
the same triangular sublattice), we have only three unit cells 
with space group P31m, P31c, and Cmc2i as shown in Fig. 
6 . Note that, the P31m cell is just a doubling of single-layer 
unit cell, and the P31c structure differ from the P31m struc¬ 
ture by the position of Ru hollow sites, so that half of Ru sites 
avoid sitting on top of Ru sites in the neighboring layer as 
shown in bottom panels of Fig. 6. Finally, the Cmc2i struc¬ 
ture differs from other unit cells by anti-cyclic stacking of ev¬ 
ery other RuCls layer as shown in the stacking sequence in the 
figure, which can be obtained by applying mirror operation to 
every other RuCls layers. 

Structure optimizations were performed including SOC, 
and Table II shows the optimized lattice constants with respect 
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Lattice constant a (in A) 


FIG. 5. (Color online) (a) Schematic figure representing the zigzag 
magnetic order in the Ru honeycomb plane. Note that, the moments 
are confined on the ac plane, where the angle of the moments with 
respect to the a-axis is denoted as 0. (b) Evolution of 0 for the two 
zigzag-ordered chains in the monoclinic unit cell as a function of a, 
obtained from OpenMX (blue) and VASP (red) codes. Structures 
with a = 5.84, 5.96, and 6.023A are Case I, II, and III, respectively, 
and angles from the structures are marked as filled symbols. Rest of 
the results are obtained from interpolation between the three struc¬ 
tures, marked as empty symbols. 


P3il2 C2/m P31m P31c Cmc2i 

Lattce constants 


a/ao 

0.984 

0.981 

0.986 

0.985 

0.984 

b/ao 

0.984 

0.986 

0.986 

0.985 

0.983 

c/co 

1.014 

1.013 

1.005 

1.007 

1.014 

AE / f.u. 

(in meV) 

1.4 

1.4 

0.0 

2.8 

2.5 

DOS at Ef 

(in states / eV / f.u.) 

9.2 

7.9 

6.0 

10.8 

8.5 


TABLE II. Optimized lattice constants, relative total energies (AE) 
per formula unit (f.u.), and densities of states (DOS) at the Eermi 
level for five stacking unit cells. Values are obtained using PBEsol 
functional and including SOC, but without electron interactions. 


to experimentally reported lattice constants ao = 5.96A and 
Co = 17.2A and their relative total energies. Note that, struc¬ 
tures without threefold symmetry — monoclinic (72/m and 
orthorhombic Cmc2i — shows slightly different a/ao and 
b/bo. Among the five different structures, the P31m structure 
yields the lowest energy. The P3il2 and (72/m structures 
are closer in energy by 1.4 meV / f.u., and for the other phases 


P3il2 (72/m P31m P31c Cmc2i 


Lattice constants 

a/ao 

1.011 

1.011 

1.010 

1.011 

1.010 

b/bo 

1.006 

1.006 

1.006 

1.006 

1.006 

c/co 

1.041 

1.043 

1.067 

1.039 

1.056 

AE / f.u. (meV) 

cFM 

0.4 

0.1 

3.7 

0.8 

0.0 

cAF 

0.4 

0.2 

4.1 

0.9 

0.4 


TABLE III. Optimized lattice constants for five stacking unit cells 
with using PBEsol functional and including SOC, Peff and mag¬ 
netism. a, 6, and c are the optimized monoclinic lattice constants 
(shown in Eig. 1) with ao, bo, and co being their experimentally 
observed values, respectively^. 


energy differences are less than 3 meV / f.u. compared to the 
the P31m structure. The lowest energy of the P31m struc¬ 
ture can be attributed to the lager kinetic energy gain originat¬ 
ing from the larger band dispersion along the c-direction com¬ 
pared to other structures. This is reflected in the lower DOS 
of the P31m cell at the Fermi level compared to other struc¬ 
tures, as shown in Table II and Fig. 7. Fig. 7 presents total 
DOS for the five structures in the presence of SOC. Compared 
to the single-layer result depicted as grey shade in the figure, 
layer stacking yields pronounced peaks near the Fermi level 
except the P31m structure in the results without SOC (not 
shown) due to the presence of flat bands along the c-direction 
at the Fermi level. Inclusion of SOC smoothes the peaks, but 
the gross feature remains the same as shown in Fig. 7, so 
resulting in higher DOS at the Fermi level except the P31m 
structure as shown in Table II. Note that, Stoner-type ferro¬ 
magnetic (FM) instability is also observed, but in this study 
we concentrate on the experimentally observed zigzag mag¬ 
netic order as discussed in the next section. 


V. STACKING WITH ZIGZAG MAGNETIC ORDER 

Now we present the stacking results that include the on¬ 
site Coulomb interaction and magnetism. Fig. 8 shows 10 
trial structural and magnetic configurations, where the direc¬ 
tion of magnetic moments in each layer is the same with the 
single-layer result in Sec. III. Fixing the in-plane zigzag or¬ 
der, we chose two interlayer magnetic configurations that we 
denote as cFM and cAF hereafter. As shown in Fig. 8, in 
the cFM configuration the zigzag-ordered layers are stacked 
along the c-direction so that the FM zigzag chains in adjacent 
layers become closer in distance, while in the cAF configura¬ 
tion the moments on one Ru layer are flipped. Note that, there 
can be additional magnetic stacking orders due to the three¬ 
fold rotational degree of freedom for each single-layer zigzag 
order, — three different direction for FM zigzag chains — 
and in this work we chose the simplest conflguration com¬ 
mensurate to the monoclinic unit cell (shown in Fig. 2(a)) 
for each structure. Structural optimizations were done first 
by varying c-axis with fixing a-lattice constants determined in 
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FIG. 6. (Color online) Five different unit cells with two- and three-layer periodicity along c-direction. Upper and lower panels show the side 
view of the unit cell and schematic top view of Ru honeycomb stacking, respectively. 


P3112 





Cmc2i 

-1.2 -0.9 -0.6 -0.3 0.0 0.3 

Energy (eV) 

FIG. 7. Densities of states (DOS) for different a-RuCU structures 
including SOC in the absence of Ueff and magnetism. Grey shade 
shows DOS of single-layer RuCU multiplied by 0.5 as a reference. 


the single-layer calculation, and later fully optimizing a and 
c axis constants and internal coordinates. Note that, symme¬ 
try constraints are lost during the full optimizations including 
the Coulomb interaction and magnetism. As a result, the opti¬ 
mized structures slightly deviate from the original space group 
symmetries, where the deviation develops in Cl positions with 
its size about 1% for each internal coordinate compared to 
the lattice constants. Note also that, structural optimization 


P3il2 C2/m 



P31m P31c Cmc2i 

~CX5 CO 00 

-00 00 00 

FIG. 8. (Color online) 10 trial magnetic configurations with in-plane 
zigzag order, where red and blue symbols depicting Ru sites with 
antiparallel magnetic moments to each other. 


for each stacking with different magnetic configuration (either 
cFM or cAF configurations in Fig. 8) yielded negligible dif¬ 
ferences. Ah of the configurations become insulator with the 
gap of ^ 1 eV between the lower and upper Hubbard bands at 
Ueff = 2 eV. DOS for the resulting phases are almost identical 
to the one from single-layer calculation^^ and show no signif¬ 
icant difference compared to each other, so we do not present 
the DOS plots here. 

Table III shows the optimization results. Compared to the 
results without stacking and magnetism, a few differences can 
be noticed; i) Energy differences between structures are less 
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FIG. 9. (Color online) Two distinct NN bonds — Z- and X/Y-bonds 
— in the presence of zigzag magnetic order, and the dominant in¬ 
terlayer hopping channels in C2 /?ti structure. Note that, although 
the largest interlayer hopping channels is depicted as green thick and 
solid arrows, magnitudes of all of the hopping channels in the figure 
are comparable to each other. 


than 1 meV per f.u. except the P31m structure, which is 
higher in energy by ^ 4.0 meV / f.u. compared to other struc¬ 
tures. Note that, the P31m structure showed the lowest en¬ 
ergy in the calculation without Peff and magnetism. With Peff 
and magnetism introduced, gap is fully opened for all of the 
structures and the relative energy gain in the P31m structure 
due to the c-axis dispersion (discussed in Sec. IV) becomes 
smaller, ii) Energy differences between cFM and cAF config¬ 
urations are smaller than 0.1 meV /f.u. for the P3i 12, C2/m, 
and P31c structures, and for the P31c and Cmc2i stackings 
the differences are about 0.4 meV / f.u.. Such small energy 
differences can be attributed to weak interlayer exchange in¬ 
teractions, which will be discussed later in the last paragraph 
of this section, iii) Fattice constants are increased by 2 to 3 
% compared to the results without Peff • iv) Small monoclinic 
distortion, which manifests itself by the difference of a/ao 
and b/bo (and negative 6 in Fig. 2), happens in every struc¬ 
tures in the presence of the in-plane zigzag magnetic order. 

Except the P31m structure which is higher in energy by ^ 
4 meV / f.u. compared to other structures, the structural en¬ 
ergy differences are smaller than I meV. This result implies 
the coexistence of different structures in experimentally syn¬ 
thesizes samples. Especially, it is natural that the P3il2 and 
C2/m structures have similar total energies; their only differ¬ 
ence is the stacking of the Ru honeycomb order, which can 
be switched to each other by the ionic hopping of Ru atoms 
within the RuCls layers. Indeed, both were reported as the 
crystal structure of a-RuCls by different groups^’^^’^^. It is 
also interesting that, the Cmc2i structure (with cFM order) 
shows the lowest energy, which can be transformed into other 
structures by applying mirror operations to every other RuCls 
layers. One can speculate that the Cmc2i structure forms in 
high temperature regime and freeze below T ^ 150 K, where 
an anomalous behavior in magnetic susceptibility observed^’^, 
so contributing to the magnetic peak with two-layer periodic¬ 
ity in poly crystalline samples below T/vi — 14 

Finally, we comment on the interlayer exchange interac¬ 
tions. Major interlayer hopping channels are shown in Fig. 
9, where the largest channel is depicted as green solid arrow 
while others are represented as dashed/dotted arrows. Note 
that, value of the largest interlayer t 2 g hopping term is about 


35 meV, and magnitudes of other channels depicted in the fig¬ 
ure are comparable to the largest one; about 20 to 30 meV. 
The interlayer exchange Heisenberg term is roughly estimated 
to be J = ^ 0.05 meV for the jeff = 1/2 pseu¬ 

dospins. This value is two orders-of-magnitude smaller than 
the previously estimated in-plane exchange interactions in a- 
RuCls^’^^ and is also consistent with the small energy differ¬ 
ences between the cFM and cAF phases discussed above. 


VI. DISCUSSION 

The relative energies between different stacking order de¬ 
pends on the electronic structures of each system in our re¬ 
sults, especially whether the system becomes fully insulating 
or not. Given that a-RuCls remains insulating in the paramag¬ 
netic phase above Tati with I eV of optical gap"^’^^, we spec¬ 
ulate that the four stacking orders — P3il2, C2/m, P31c, 
and Cmc2i — are almost degenerate as discussed in Sec. V. 

The change of hopping integrals and exchange interaction 
terms after the structure optimization show that the physics of 
a-RuCls is sensitive the NN Ru-Ru distance and Cl position. 
For example, the strength of the Kitaev and F terms are sig¬ 
nificantly modified by the Ru-Ru and Ru-Cl distances. This 
implies that, even a small amount of epitaxial tensile strain by 
1 % or uniaxial pressure perpendicular to the layer can signif¬ 
icantly enhance the Kitaev term and push the system closer 
to the Kitaev limit. On the other hand, hydrostatic pressure 
or compressive strain can increase the ts term by decreasing 
the Ru-Ru distance. This reduces the FM Kitaev term and 
drive the effective model to the highly frustrated F-dominated 
regime. In addition, presence of the negative F' term due to 
the trigonal distortion can stablize the zigzag-ordered phase 
as discussed in previous study^^. Effects of the monoclnic 
bond disproportionation'^'^ is another factor that can change 
the magnetism. In this regard, full experimental structure de¬ 
termination including precise atomic positions and stacking 
order would be important for future studies. 

In summary, structural properties of a-RuCls from ab- 
initio calculations are presented in this study. SOC is found 
to prevent the Ru dimerization in the Ru honeycomb layers, 
and the presence of in-plane zigzag magnetic order further 
gives small monoclinic distortion. The relation between the 
hopping integrals and exchange interactions to the structure 
is also discussed. Total energy comparison between different 
RuCls stacking orders yields the Cmc2i and C2/m struc¬ 
tures to be the almost degenerate ground state structures, and 
P3il2 structure to be comparable in energy; energy differ¬ 
ences smaller than 0.4 meV per formula unit. In-plane ex¬ 
change interactions are found to be sensitive to the structural 
distortions, and the jeff = 1/2 pseudospin model is domi¬ 
nated by the FM Kitaev terms in the optimized structures with 
the presence of Peff, similar to the two- and three-dimensional 
honeycomb iridates'^^"'^^ . As expected, interlayer exchange in¬ 
teractions are estimated to be weak compared to the in-plane 
exchange interactions, so this system can be a good platform 
in studying frustrated two-dimensional magnetism. 

Note added — After the completion of the manuscript, we 
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became aware of the experimental work by Johnson and co- 
workers'^^, which reports monoclinic C2lm crystal structure 
and the in-plane zigzag magnetic configuration with antifer¬ 
romagnetic interplanar order below T/v ^ 13 K. 
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Appendix A: van der Waals calculation 

In this Appendix we compare the results with using dif¬ 
ferent exchange-correlation functionals including vdW in¬ 
teractions. Four functionals are considered; PBE, PBEsol, 
vdW-DF24® and vdW-optB86b^o, where vdW-DF2 and vdW- 
optB86b functionals showed accuracies comparable to RPA 


calculations in layered and bulk systems, respectively. Here 
SOC, Ueff, and magnetism are not included. 

Fig. A1 shows relative energies versus c-lattice constant 
with fixed a = ao for the results with four functionals, where 
C2/m stacking order is not considered. Except PBE, which 
yields unreasonably large value of c, other three functionals 
yields P31m and P31c as configurations with the lowest and 
second lowest energy. Compared to PBEsol, vdW functionals 
tend to yield steeper energy curve away from the optimum c 
valuex and higher energy for P3il2 phase. 

Table A1 shows the results from full lattice optimizations. 
Except the change of a-lattice constants, where vdW-DF2 
results yields 3% enhancement of a value, the features are 
qualitatively similar to the results in Fig. Al. P31m is 
still the most favored configuration, and optimized c-lattice 
constants do not change significantly from the values in Fig. 
Al. It is notable that the vdW results give high energies for 
P3il2 and Cmc2i phases, which were the favored phases in 
PBEsoUSOC-FPeff calculations. 

Compared to the vdW functionals, PBEsol yields reason¬ 
able estimates of total energy and lattice constants, although 
quantitative differences can be noticed. Since test calculations 
on combining vdW functionals and DFT-fSOC-hP, which is 
crucial in understanding physics of RuCls, have not been done 
yet, in this study PBEsol functional is employed for the rest 
of the calculations. 
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FIG. Al. (Color online) Total energy vs. c-lattice constant plots with fixed a = ao = 5.96A for different crystal structure and exchange- 
correlation functionals. From left to right, results with PBE, PBEsol, vdW-DE2 and vdW-optB86b functionals are shown. 
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